Introduction

This is a reanalysis of the NC13955 population QTL analysis. This is done to produce a single, specific QTL analysis for publication.

#set working directory
setwd("C:/Users/zwinn/OneDrive/Publication/Manuscripts In Progress/FHB QTL/Final Formatted Data and Analysis")

#read in data
pheno<-read.csv("NC13955 Combined Data Final.csv")

#make sure everything is the right value
pheno[,1:7]<-lapply(pheno[,1:7], as.factor)
pheno[,8:ncol(pheno)]<-lapply(pheno[,8:ncol(pheno)], as.numeric)

#read in genetic information
geno<-read.vcf("NC13-20076xGA06493-13LE6_filt.vcf.gz", convert.chr = F)
## ped stats and snps stats have been set. 
## 'p' has been set. 
## 'mu' and 'sigma' have been set.
#convert genotype matrix to format for rqtl 
geno<-format_qtlmap_geno(geno, 
                         par_a="13955-GA06493-13LE6",
                         par_b="13955-NC13-20076",
                         rm_het=T, 
                         rm_miss=T, 
                         include_pars=T, 
                         out_fmt="rqtl")$abh

#replace heterozygous calls with missing data
geno<-replace(geno, geno=="H", "-")

Calculate BLUEs within and across environments

#set up lists for loop
env<-levels(pheno$ENV)
traits<-colnames(pheno)[8:12]

#create dataframe for means out
BLUEs_loc<-data.frame(GENOTYPE=unique(pheno$GENOTYPE))

#check normality
for (i in traits){
  hist(pheno[[i]], main=i)
}

#run models for each environment
for (j in env){
  for (i in traits){
  location<-pheno[pheno$ENV==j,] #filter out each environment
  print(paste("------- Analyzing trait", i, "in", j,"-------"))
  if((sum(is.na(location[[i]]==T)))>200){  #check if there is data
    print("Data not found, moving to the next trait")
    }else
      if (i=="FDK"){
        
        print(paste("Analyzing with GMLM with Inverse Gaussian distribution"))
        
        #Run Model
        mlm<-asreml(fixed = location[,i]~1+GENOTYPE, 
              random = ~ (REP), 
              residual = ~ idv(units), 
              data = location, 
              family= asr_inverse.gaussian())
        
        #Check significance
        print(wald(mlm))
        
        #print resid
        qqnorm(mlm$residuals)
        
        #Pull BLUEs  
        beta_hat<-mlm$coefficients$fixed
        beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
        beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
        BLUEs<-as.data.frame(beta_0+beta_1)
        rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
        colnames(BLUEs)<-paste(i, "_", j, sep = "")
        BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
        BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
        
        #Remove Junk
        remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
        
      }else
        if(i=="DON"){
        
          print(paste("Analyzing with GMLM with Inverse Gaussian distribution"))
        
          #Run Model
          mlm<-asreml(fixed = location[,i]~1+GENOTYPE,
                      random = ~ (REP), 
                      residual = ~ idv(units), 
                      data = location, 
                      family= asr_inverse.gaussian())
         
          #Check significance
          print(wald(mlm))
          
          #print resid
          qqnorm(mlm$residuals)
          
          #Pull BLUEs  
          beta_hat<-mlm$coefficients$fixed
          beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
          beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
          BLUEs<-as.data.frame(beta_0+beta_1)
          rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
          colnames(BLUEs)<-paste(i, "_", j, sep = "")
          BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
          BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
          
          #Remove Junk
          remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
        
          
        }else{
          
          print(paste("Analyzing with MLM with Gaussian distribution"))
        
          #Run Model
          mlm<-asreml(fixed = location[,i]~1+GENOTYPE, 
                      random = ~ (REP), 
                      residual = ~ idv(units), 
                      data = location)
         
          #Check significance
          print(wald(mlm))
          
          #print resid
          qqnorm(mlm$residuals)
          
          #Pull BLUEs  
          beta_hat<-mlm$coefficients$fixed
          beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
          beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
          BLUEs<-as.data.frame(beta_0+beta_1)
          rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
          colnames(BLUEs)<-paste(i, "_", j, sep = "")
          BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
          BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
          
          #Remove Junk
          remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
        }
  }
}
## [1] "------- Analyzing trait HD in KIN2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FHB in KIN2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:08 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -228.014           1.0    149 13:35:08    0.0 (1 restrained)
##  2      -181.290           1.0    149 13:35:08    0.0 (1 restrained)
##  3      -133.124           1.0    149 13:35:09    0.0 (1 restrained)
##  4      -104.473           1.0    149 13:35:09    0.0 (1 restrained)
##  5       -93.401           1.0    149 13:35:09    0.0 (1 restrained)
##  6       -92.100           1.0    149 13:35:09    0.0
##  7       -92.070           1.0    149 13:35:09    0.0
##  8       -92.070           1.0    149 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    7676.9         7676.9 < 2.2e-16 ***
## GENOTYPE      184    1420.4         1420.4 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FDK in KIN2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -850.740           1.0    158 13:35:09    0.0
##  2      -757.507           1.0    158 13:35:09    0.0
##  3      -656.700           1.0    158 13:35:09    0.0
##  4      -590.195           1.0    158 13:35:09    0.0
##  5      -558.335           1.0    158 13:35:09    0.0
##  6      -551.609           1.0    158 13:35:09    0.0
##  7      -551.069           1.0    158 13:35:09    0.0
##  8      -551.063           1.0    158 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1   2288.37        2288.37 < 2.2e-16 ***
## GENOTYPE      184    862.08         862.08 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in KIN2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -819.956           1.0    158 13:35:09    0.0
##  2      -737.235           1.0    158 13:35:09    0.0
##  3      -648.235           1.0    158 13:35:09    0.0
##  4      -590.321           1.0    158 13:35:09    0.0
##  5      -563.424           1.0    158 13:35:09    0.0
##  6      -558.168           1.0    158 13:35:09    0.0
##  7      -557.780           1.0    158 13:35:09    0.0
##  8      -557.771           1.0    158 13:35:09    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    233.26         233.26 < 2.2e-16 ***
## GENOTYPE      184    967.91         967.91 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in KIN2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -641.368           1.0    232 13:35:09    0.0 (1 restrained)
##  2      -537.434           1.0    232 13:35:09    0.0 (1 restrained)
##  3      -427.349           1.0    232 13:35:09    0.0 (1 restrained)
##  4      -357.531           1.0    232 13:35:09    0.0 (1 restrained)
##  5      -326.764           1.0    232 13:35:09    0.0 (1 restrained)
##  6      -321.594           1.0    232 13:35:09    0.0
##  7      -321.337           1.0    232 13:35:09    0.0
##  8      -321.336           1.0    232 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1   1021841        1021841 < 2.2e-16 ***
## GENOTYPE      192      1327           1327 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FHB in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -882.862           1.0    231 13:35:09    0.0
##  2      -693.549           1.0    231 13:35:09    0.0
##  3      -486.358           1.0    231 13:35:09    0.0
##  4      -344.914           1.0    231 13:35:09    0.0
##  5      -271.571           1.0    231 13:35:09    0.0
##  6      -252.398           1.0    231 13:35:09    0.0
##  7      -249.968           1.0    231 13:35:09    0.0
##  8      -249.902           1.0    231 13:35:09    0.0
##  9      -249.902           1.0    231 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    352.28         352.28 < 2.2e-16 ***
## GENOTYPE      193    732.77         732.77 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FDK in KIN2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1169.203           1.0    233 13:35:09    0.0
##  2     -1043.001           1.0    233 13:35:09    0.0
##  3      -906.989           1.0    233 13:35:09    0.0
##  4      -818.037           1.0    233 13:35:09    0.0
##  5      -776.170           1.0    233 13:35:09    0.0
##  6      -767.534           1.0    233 13:35:09    0.0
##  7      -766.662           1.0    233 13:35:09    0.0
##  8      -766.560           1.0    233 13:35:09    0.0
##  9      -766.546           1.0    233 13:35:09    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    541.64         541.64 < 2.2e-16 ***
## GENOTYPE      193   1141.84        1141.84 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in KIN2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -671.414           1.0    182 13:35:09    0.0
##  2      -600.399           1.0    182 13:35:09    0.0
##  3      -525.043           1.0    182 13:35:09    0.0
##  4      -477.866           1.0    182 13:35:09    0.0
##  5      -457.810           1.0    182 13:35:09    0.0
##  6      -454.704           1.0    182 13:35:09    0.0
##  7      -454.540           1.0    182 13:35:09    0.0
##  8      -454.534           1.0    182 13:35:09    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    190.79         190.79 < 2.2e-16 ***
## GENOTYPE      193   1325.51        1325.51 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -860.120           1.0    234 13:35:09    0.0 (1 restrained)
##  2      -743.455           1.0    234 13:35:09    0.0
##  3      -619.080           1.0    234 13:35:09    0.0
##  4      -538.702           1.0    234 13:35:09    0.0
##  5      -501.992           1.0    234 13:35:09    0.0
##  6      -495.193           1.0    234 13:35:09    0.0
##  7      -494.779           1.0    234 13:35:09    0.0
##  8      -494.777           1.0    234 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1     86974          86974 < 2.2e-16 ***
## GENOTYPE      193      1211           1211 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait HD in LW2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -691.787           1.0    192 13:35:09    0.0 (1 restrained)
##  2      -559.886           1.0    192 13:35:09    0.0
##  3      -416.268           1.0    192 13:35:09    0.0
##  4      -319.913           1.0    192 13:35:09    0.0
##  5      -271.918           1.0    192 13:35:09    0.0
##  6      -260.652           1.0    192 13:35:09    0.0
##  7      -259.517           1.0    192 13:35:09    0.0
##  8      -259.498           1.0    192 13:35:09    0.0
##  9      -259.498           1.0    192 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    300979         300979 < 2.2e-16 ***
## GENOTYPE      186       804            804 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FHB in LW2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -317.568           1.0    186 13:35:09    0.0 (1 restrained)
##  2      -247.435           1.0    186 13:35:09    0.0 (1 restrained)
##  3      -173.866           1.0    186 13:35:09    0.0 (1 restrained)
##  4      -128.702           1.0    186 13:35:09    0.0 (1 restrained)
##  5      -110.014           1.0    186 13:35:09    0.0 (1 restrained)
##  6      -107.354           1.0    186 13:35:09    0.0
##  7      -107.261           1.0    186 13:35:09    0.0
##  8      -107.261           1.0    186 13:35:09    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    9558.3         9558.3 < 2.2e-16 ***
## GENOTYPE      185    1358.2         1358.2 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FDK in LW2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:09 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1337.441           1.0    197 13:35:10    0.0
##  2     -1153.277           1.0    197 13:35:10    0.0
##  3      -950.889           1.0    197 13:35:10    0.0
##  4      -811.154           1.0    197 13:35:10    0.0
##  5      -736.793           1.0    197 13:35:10    0.0
##  6      -715.957           1.0    197 13:35:10    0.0
##  7      -712.894           1.0    197 13:35:10    0.0
##  8      -712.779           1.0    197 13:35:10    0.0
##  9      -712.779           1.0    197 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1   1430.25        1430.25 < 2.2e-16 ***
## GENOTYPE      186    600.11         600.11 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in LW2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1414.045           1.0    197 13:35:10    0.0
##  2     -1216.349           1.0    197 13:35:10    0.0
##  3      -998.602           1.0    197 13:35:10    0.0
##  4      -847.293           1.0    197 13:35:10    0.0
##  5      -765.488           1.0    197 13:35:10    0.0
##  6      -741.388           1.0    197 13:35:10    0.0
##  7      -737.166           1.0    197 13:35:10    0.0
##  8      -736.727           1.0    197 13:35:10    0.0
##  9      -736.649           1.0    197 13:35:10    0.0
## 10      -736.640           1.0    197 13:35:10    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1     61.23          61.23 5.107e-15 ***
## GENOTYPE      186    569.20         569.20 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in LW2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -972.328           1.0    222 13:35:10    0.0
##  2      -772.024           1.0    222 13:35:10    0.0
##  3      -552.110           1.0    222 13:35:10    0.0
##  4      -400.678           1.0    222 13:35:10    0.0
##  5      -320.581           1.0    222 13:35:10    0.0
##  6      -298.500           1.0    222 13:35:10    0.0
##  7      -295.365           1.0    222 13:35:10    0.0
##  8      -295.256           1.0    222 13:35:10    0.0
##  9      -295.256           1.0    222 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1     40395          40395 < 2.2e-16 ***
## GENOTYPE      192       635            635 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FHB in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -482.340           1.0    223 13:35:10    0.0 (1 restrained)
##  2      -381.960           1.0    223 13:35:10    0.0 (1 restrained)
##  3      -275.882           1.0    223 13:35:10    0.0
##  4      -208.294           1.0    223 13:35:10    0.0
##  5      -178.416           1.0    223 13:35:10    0.0 (1 restrained)
##  6      -173.362           1.0    223 13:35:10    0.0 (1 restrained)
##  7      -173.108           1.0    223 13:35:10    0.0
##  8      -173.107           1.0    223 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    6294.6         6294.6 < 2.2e-16 ***
## GENOTYPE      192    1294.8         1294.8 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FDK in LW2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1079.621           1.0    218 13:35:10    0.0
##  2      -958.974           1.0    218 13:35:10    0.0
##  3      -828.869           1.0    218 13:35:10    0.0
##  4      -743.660           1.0    218 13:35:10    0.0
##  5      -703.515           1.0    218 13:35:10    0.0
##  6      -695.396           1.0    218 13:35:10    0.0
##  7      -694.784           1.0    218 13:35:10    0.0
##  8      -694.776           1.0    218 13:35:10    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    444.85         444.85 < 2.2e-16 ***
## GENOTYPE      192   1070.34        1070.34 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait DON in LW2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -356.315           1.0    106 13:35:10    0.0
##  2      -302.485           1.0    106 13:35:10    0.0
##  3      -244.630           1.0    106 13:35:10    0.0
##  4      -207.098           1.0    106 13:35:10    0.0
##  5      -189.805           1.0    106 13:35:10    0.0
##  6      -186.524           1.0    106 13:35:10    0.0
##  7      -186.313           1.0    106 13:35:10    0.0
##  8      -186.312           1.0    106 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    438.52         438.52 < 2.2e-16 ***
## GENOTYPE      184    851.54         851.54 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -902.787           1.0    223 13:35:10    0.0 (1 restrained)
##  2      -765.630           1.0    223 13:35:10    0.0
##  3      -618.001           1.0    223 13:35:10    0.0
##  4      -520.171           1.0    223 13:35:10    0.0
##  5      -472.821           1.0    223 13:35:10    0.0
##  6      -462.545           1.0    223 13:35:10    0.0
##  7      -461.669           1.0    223 13:35:10    0.0
##  8      -461.658           1.0    223 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    130425         130425 < 2.2e-16 ***
## GENOTYPE      191       976            976 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait HD in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FHB in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FDK in VIR2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1020.821           1.0    195 13:35:10    0.0
##  2      -922.179           1.0    195 13:35:10    0.0
##  3      -816.221           1.0    195 13:35:10    0.0
##  4      -747.594           1.0    195 13:35:10    0.0
##  5      -716.091           1.0    195 13:35:10    0.0
##  6      -710.172           1.0    195 13:35:10    0.0
##  7      -709.799           1.0    195 13:35:10    0.0
##  8      -709.796           1.0    195 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    1028.7         1028.7 < 2.2e-16 ***
## GENOTYPE      186    1080.5         1080.5 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in VIR2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -576.052           1.0    196 13:35:10    0.0
##  2      -535.425           1.0    196 13:35:10    0.0
##  3      -493.851           1.0    196 13:35:10    0.0
##  4      -470.318           1.0    196 13:35:10    0.0
##  5      -462.305           1.0    196 13:35:10    0.0
##  6      -461.553           1.0    196 13:35:10    0.0
##  7      -461.504           1.0    196 13:35:10    0.0
##  8      -461.500           1.0    196 13:35:10    0.0
## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.

## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    148.11         148.11 < 2.2e-16 ***
## GENOTYPE      186   2082.22        2082.22 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in VIR2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:10 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -611.866           1.0    234 13:35:10    0.0 (1 restrained)
##  2      -495.363           1.0    234 13:35:10    0.0
##  3      -371.420           1.0    234 13:35:10    0.0
##  4      -291.397           1.0    234 13:35:10    0.0
##  5      -254.889           1.0    234 13:35:10    0.0
##  6      -248.148           1.0    234 13:35:10    0.0
##  7      -247.741           1.0    234 13:35:10    0.0
##  8      -247.738           1.0    234 13:35:10    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1   2559200        2559200 < 2.2e-16 ***
## GENOTYPE      194      1219           1219 < 2.2e-16 ***
## residual (MS)             1                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FHB in VIR2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:11 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -614.597           1.0    233 13:35:11    0.0 (1 restrained)
##  2      -492.704           1.0    233 13:35:11    0.0
##  3      -362.606           1.0    233 13:35:11    0.0
##  4      -278.034           1.0    233 13:35:11    0.0
##  5      -238.879           1.0    233 13:35:11    0.0
##  6      -231.352           1.0    233 13:35:11    0.0
##  7      -230.856           1.0    233 13:35:11    0.0
##  8      -230.853           1.0    233 13:35:11    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    1047.2         1047.2 < 2.2e-16 ***
## GENOTYPE      194    1163.7         1163.7 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait FDK in VIR2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:11 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1116.367           1.0    228 13:35:11    0.0
##  2      -989.165           1.0    228 13:35:11    0.0
##  3      -851.930           1.0    228 13:35:11    0.0
##  4      -761.943           1.0    228 13:35:11    0.0 (1 restrained)
##  5      -719.476           1.0    228 13:35:11    0.0 (1 restrained)
##  6      -710.884           1.0    228 13:35:11    0.0 (1 restrained)
##  7      -710.252           1.0    228 13:35:11    0.0 (1 restrained)
##  8      -710.246           1.0    228 13:35:11    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    2234.4         2234.4 < 2.2e-16 ***
## GENOTYPE      194    1085.2         1085.2 < 2.2e-16 ***
## residual (MS)           1.0                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait DON in VIR2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:11 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1      -327.582           1.0    225 13:35:11    0.0 (1 restrained)
##  2      -272.817           1.0    225 13:35:11    0.0
##  3      -215.680           1.0    225 13:35:11    0.0
##  4      -182.480           1.0    225 13:35:11    0.0
##  5      -170.610           1.0    225 13:35:11    0.0
##  6      -169.501           1.0    225 13:35:11    0.0
##  7      -169.485           1.0    225 13:35:11    0.0
##  8      -169.485           1.0    225 13:35:11    0.0
## Wald tests for fixed effects.
## Response: location[, i]
## Terms added sequentially; adjusted for those above.
## 
##                Df Sum of Sq Wald statistic Pr(Chisq)    
## (Intercept)     1    552.94         552.94 < 2.2e-16 ***
## GENOTYPE      194   2020.25        2020.25 < 2.2e-16 ***
## residual (MS)          1.00                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "------- Analyzing trait PH in VIR2020 -------"
## [1] "Data not found, moving to the next trait"

Multi-Environmental Models

#define traits
traits<-colnames(pheno)[8:12]

#make dataframe for predictions
BLUEs_mem<-data.frame(GENOTYPE=unique(pheno$GENOTYPE))

#run MEMLM loop
for (i in traits){
  print(paste("------- Analyzing trait", i, "-------"))
  
  siteyears<-pheno %>% 
    select(ENV, YEAR, LOC, all_of(i)) %>% 
    drop_na() %>% 
    distinct(YEAR)
  
  siteyears<-as.numeric(count(siteyears))
  
  if(((siteyears)>=2)){
    
    print(paste(i, "has multiple years of data"))
    
    if(i=="FDK"){
      #run models    
      mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
                  random = ~ ENV+
                    GENOTYPE:ENV,
                  residual = ~idv(units),
                  data = pheno,
                  maxit=75,
                  family= asr_inverse.gaussian())
      
      #print resid
      qqnorm(mlm$residuals)
      
      #Pull BLUEs  
      beta_hat<-mlm$coefficients$fixed
      beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
      beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
      BLUEs<-as.data.frame(beta_0+beta_1)
      rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
      colnames(BLUEs)<-paste(i, "_ME", sep = "")
      BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
      BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
      
      #Remove Junk
      remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
    }else
      if(i=="DON"){
        #run models    
        mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
                    random = ~ ENV+
                      GENOTYPE:ENV,
                    data = pheno,
                    maxit=75,
                    family= asr_inverse.gaussian())
        
        #print resid
        qqnorm(mlm$residuals)
        
        #Pull BLUEs  
        beta_hat<-mlm$coefficients$fixed
        beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
        beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
        BLUEs<-as.data.frame(beta_0+beta_1)
        rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
        colnames(BLUEs)<-paste(i, "_ME", sep = "")
        BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
        BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
        
        #Remove Junk
        remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
      }else{
        #run models    
        mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
                    random = ~ENV+
                      GENOTYPE:ENV,
                    residual = ~idv(units),
                    data = pheno,
                    maxit=75)
        
        #print resid
        qqnorm(mlm$residuals)
        
        #Pull BLUEs  
        beta_hat<-mlm$coefficients$fixed
        beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
        beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
        BLUEs<-as.data.frame(beta_0+beta_1)
        rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
        colnames(BLUEs)<-paste(i, "_ME", sep = "")
        BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
        BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
        
        #Remove Junk
        remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
      }
  }else{
    
    print(paste(i, "has a single year of data"))
    
    #run model
    mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
                random = ~ENV+
                  GENOTYPE:ENV,
                residual = ~units,
                data = pheno,
                maxit=75)
    
    #print resid
    qqnorm(mlm$residuals)
    
    #Pull BLUEs  
    beta_hat<-mlm$coefficients$fixed
    beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
    beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
    BLUEs<-as.data.frame(beta_0+beta_1)
    rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
    colnames(BLUEs)<-paste(i, "_ME", sep = "")
    BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
    BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
    
    #Remove Junk
    remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
    }
}
## [1] "------- Analyzing trait HD -------"
## [1] "HD has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:12 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -2381.807           1.0   1451 13:35:12    0.0 (1 restrained)
## Log-likelihood decreased to  -2678.06 ; trying with reduced updates (0.81626)
##  2     -2325.628           1.0   1451 13:35:12    0.0
##  3     -2030.370           1.0   1451 13:35:12    0.0
##  4     -1798.807           1.0   1451 13:35:12    0.0
##  5     -1771.640           1.0   1451 13:35:12    0.0
##  6     -1763.747           1.0   1451 13:35:12    0.0
##  7     -1760.076           1.0   1451 13:35:12    0.0
##  8     -1758.663           1.0   1451 13:35:12    0.0
##  9     -1758.233           1.0   1451 13:35:12    0.0
## 10     -1758.157           1.0   1451 13:35:12    0.0
## 11     -1758.153           1.0   1451 13:35:12    0.0

## [1] "------- Analyzing trait FHB -------"
## [1] "FHB has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:12 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -3215.668           1.0   1779 13:35:12    0.0
##  2     -2618.147           1.0   1779 13:35:12    0.0
##  3     -1986.339           1.0   1779 13:35:12    0.0
##  4     -1593.414           1.0   1779 13:35:12    0.0
##  5     -1426.382           1.0   1779 13:35:12    0.0
##  6     -1399.236           1.0   1779 13:35:12    0.0
##  7     -1397.875           1.0   1779 13:35:12    0.0
##  8     -1397.871           1.0   1779 13:35:12    0.0

## [1] "------- Analyzing trait FDK -------"
## [1] "FDK has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:12 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -12793.71           1.0   2174 13:35:12    0.0
##  2     -11133.10           1.0   2174 13:35:12    0.0
##  3      -9316.39           1.0   2174 13:35:12    0.0
##  4      -8075.14           1.0   2174 13:35:12    0.0
##  5      -7421.83           1.0   2174 13:35:12    0.0
##  6      -7226.13           1.0   2174 13:35:12    0.0
##  7      -7172.62           1.0   2174 13:35:12    0.0
##  8      -7152.58           1.0   2174 13:35:12    0.0
##  9      -7143.52           1.0   2174 13:35:12    0.0
## 10      -7139.69           1.0   2174 13:35:12    0.0
## 11      -7138.33           1.0   2174 13:35:12    0.0
## 12      -7138.00           1.0   2174 13:35:12    0.0
## 13      -7137.97           1.0   2174 13:35:12    0.0
## 14      -7137.97           1.0   2174 13:35:12    0.0

## [1] "------- Analyzing trait DON -------"
## [1] "DON has multiple years of data"
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:13 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -6247.703       137.370   2001 13:35:13    0.0
##  2     -6221.203       126.069   2001 13:35:13    0.0
##  3     -6199.139       112.958   2001 13:35:13    0.0
##  4     -6192.242       105.175   2001 13:35:13    0.0
##  5     -6191.155       102.685   2001 13:35:13    0.0
##  6     -6191.001       102.813   2001 13:35:13    0.0
##  7     -6190.992       102.769   2001 13:35:13    0.0

## [1] "------- Analyzing trait PH -------"
## [1] "PH has a single year of data"
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:13 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -1340.569       14.0952    649 13:35:13    0.0
##  2     -1339.579       13.7940    649 13:35:13    0.0
##  3     -1338.798       13.4275    649 13:35:13    0.0
##  4     -1338.562       13.1236    649 13:35:13    0.0
##  5     -1338.561       13.1084    649 13:35:13    0.0

BLUEs_mem<-BLUEs_mem %>%
  select(GENOTYPE, FHB_ME, FDK_ME, DON_ME, HD_ME, PH_ME)

pairs.panels(BLUEs_mem[,2:ncol(BLUEs_mem)], stars = T, lm=T)

Heritabilities

#make an empty data frame for the heritability estimates
h2s<-data.frame(Trait=character(), Type=character(), 
                Estimation=numeric(), SE=numeric())

#calculate the harmonic means for each trait
e_traits<-list()
r_traits<-list()

for(j in traits){
  a<-pheno %>% 
    drop_na(all_of(j))
  
  obs_per_ind<-aggregate(a[,j]~GENOTYPE, 
                         data=a,
                         length)
  r<-(1/mean(1/obs_per_ind$`a[, j]`))
  r_traits[j]<-r
  
  e<-c()
  for(i in unique(pheno$GENOTYPE)){
  a<-pheno %>%
    drop_na(j) %>%
    filter(GENOTYPE==all_of(i))
  a<-length(unique(a$ENV))
  a<-ifelse(a==0,0,1/a)
  e<-rbind(e,a)
  }

  e_traits[j]<-(1/mean(e)) 
}
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(j)
## 
##   # Now:
##   data %>% select(all_of(j))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `GENOTYPE == all_of(i)`.
## Caused by warning:
## ! Using `all_of()` outside of a selecting function was deprecated in tidyselect
##   1.2.0.
## ℹ See details at
##   <https://tidyselect.r-lib.org/reference/faq-selection-context.html>
for (i in traits){
  print(paste("------- Analyzing trait", i, "-------"))
  
  n_loc<-pheno %>% 
    select(ENV, YEAR, LOC, all_of(i)) %>% 
    drop_na() %>% 
    distinct(ENV)
  
  n_loc=as.numeric(nrow(n_loc))
  
  n_rep=n_loc*2
  
  n_year=pheno %>% 
    select(ENV, YEAR, LOC, all_of(i)) %>% 
    drop_na() %>% 
    distinct(YEAR)
  
  n_year<-as.numeric(count(n_year))
  
  e<-as.numeric(e_traits[i])
  r<-as.numeric(r_traits[i])
  
  if(((n_year)>=2)){
    
    print(paste(i, "has multiple years of data"))
    
    if(i=="FDK"){
      #Run Model
      mlm<-asreml(fixed = pheno[,i]~1,
                random = ~ GENOTYPE+
                  ENV+
                  GENOTYPE:ENV,
                residual = ~idv(units),
                data = pheno,
                maxit=75,
                family= asr_inverse.gaussian())
      
      #print summary of variance comps
      print(summary(mlm)$varcomp)
      
      #pull predictions
      pph2<-vpredict(mlm, 
                     h2~V2/(V2+V3+V4))
      
      emh2<-vpredict(mlm,
                     h2~V2/(V2+(V3/e)+(V4/(e*r))))
      
      #bind predictions
      a<-rbind(pph2,emh2)
      a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
      h2s<-rbind(h2s,a)
      
    }else
      if(i=="DON"){
  
        #Run Model
        mlm<-asreml(fixed = pheno[,i]~1,
                  random = ~ (GENOTYPE)+
                    ENV+
                    GENOTYPE:ENV,
                  residual = ~idv(units),
                  data = pheno,
                  maxit=75,
                  family= asr_inverse.gaussian())
        
        #print summary of variance comps
        print(summary(mlm)$varcomp)
        
        #pull predictions
        pph2<-vpredict(mlm, 
                       h2~V2/(V2+V3+V4))
        
        emh2<-vpredict(mlm,
                       h2~V2/(V2+(V3/e)+(V4/(e*r))))
        
        #bind predictions
        a<-rbind(pph2,emh2)
        a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
        h2s<-rbind(h2s,a)
        
      }else{
      
        #Run Model
        mlm<-asreml(fixed = pheno[,i]~1,
                  random = ~ (GENOTYPE)+
                    ENV+
                    GENOTYPE:ENV,
                  residual = ~idv(units),
                  data = pheno,
                  maxit=75)
        
        #print summary of variance comps
        print(summary(mlm)$varcomp)
        
        #pull predictions
        pph2<-vpredict(mlm, 
                       h2~V2/(V2+V3+V4))
        
        emh2<-vpredict(mlm,
                       h2~V2/(V2+(V3/e)+(V4/(e*r))))
        
        #bind predictions
        a<-rbind(pph2,emh2)
        a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
        h2s<-rbind(h2s,a)
        
    }
  }else{
    
    print(paste(i, "has a single year of data"))
    
    mlm<-asreml(fixed = pheno[,i]~1,
              random = ~ (GENOTYPE)+
                ENV+
                GENOTYPE:ENV,
              residual = ~idv(units),
              data = pheno,
              maxit=75)
        
    print(summary(mlm)$varcomp)
    
    pph2<-vpredict(mlm, 
                   h2~V2/(V2+V3+V4))
    
    emh2<-vpredict(mlm,
                   h2~V2/(V2+(V3/e)+(V4/(e*r))))
        
    a<-rbind(pph2,emh2)
    a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
    h2s<-rbind(h2s,a)
    }
}
## [1] "------- Analyzing trait HD -------"
## [1] "HD has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:25 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -2670.166           1.0   1647 13:35:26    0.0 (1 restrained)
## Log-likelihood decreased to  -2954.94 ; trying with reduced updates (0.89959)
##  2     -2607.889           1.0   1647 13:35:26    0.0
##  3     -2299.376           1.0   1647 13:35:26    0.0
##  4     -2040.728           1.0   1647 13:35:26    0.0
##  5     -2010.474           1.0   1647 13:35:26    0.0
##  6     -2002.298           1.0   1647 13:35:26    0.0
##  7     -1998.614           1.0   1647 13:35:26    0.0
##  8     -1997.200           1.0   1647 13:35:26    0.0
##  9     -1996.769           1.0   1647 13:35:26    0.0
## 10     -1996.694           1.0   1647 13:35:26    0.0
## 11     -1996.689           1.0   1647 13:35:26    0.0
##               component   std.error   z.ratio bound %ch
## ENV          245.739732 199.1495133  1.233946     P 0.5
## GENOTYPE       3.605674   0.4143780  8.701411     P 0.0
## GENOTYPE:ENV   0.607059   0.1179896  5.145023     P 0.0
## units!units    2.626815   0.1212908 21.657163     P 0.0
## units!R        1.000000          NA        NA     F 0.0
## [1] "------- Analyzing trait FHB -------"
## [1] "FHB has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:26 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -3876.512           1.0   1974 13:35:26    0.0
##  2     -3132.065           1.0   1974 13:35:26    0.0
##  3     -2339.901           1.0   1974 13:35:26    0.0
##  4     -1838.373           1.0   1974 13:35:26    0.0
##  5     -1616.057           1.0   1974 13:35:26    0.0
##  6     -1575.165           1.0   1974 13:35:26    0.0
##  7     -1572.247           1.0   1974 13:35:26    0.0
##  8     -1572.210           1.0   1974 13:35:26    0.0
##  9     -1572.210           1.0   1974 13:35:26    0.0
##              component  std.error   z.ratio bound %ch
## ENV          0.5111042 0.36461208  1.401775     P   0
## GENOTYPE     1.9613936 0.21721428  9.029763     P   0
## GENOTYPE:ENV 0.3003627 0.05016457  5.987546     P   0
## units!units  1.1396337 0.04955938 22.995319     P   0
## units!R      1.0000000         NA        NA     F   0
## [1] "------- Analyzing trait FDK -------"
## [1] "FDK has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:26 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -25196.97           1.0   2369 13:35:26    0.0
##  2     -20635.67           1.0   2369 13:35:26    0.0
##  3     -15530.46           1.0   2369 13:35:26    0.0 (1 restrained)
##  4     -11815.59           1.0   2369 13:35:26    0.0 (1 restrained)
##  5      -9558.00           1.0   2369 13:35:26    0.0 (1 restrained)
##  6      -8590.93           1.0   2369 13:35:26    0.0 (1 restrained)
##  7      -8136.15           1.0   2369 13:35:26    0.0 (1 restrained)
##  8      -7924.54           1.0   2369 13:35:26    0.0 (1 restrained)
##  9      -7835.23           1.0   2369 13:35:26    0.0
## 10      -7800.37           1.0   2369 13:35:26    0.0
## 11      -7795.71           1.0   2369 13:35:26    0.0
## 12      -7795.29           1.0   2369 13:35:26    0.0
## 13      -7795.25           1.0   2369 13:35:26    0.0
## 14      -7795.25           1.0   2369 13:35:26    0.0
##              component std.error   z.ratio bound %ch
## ENV           84.30977 53.697038  1.570101     P   0
## GENOTYPE     282.11502 30.766069  9.169681     P   0
## GENOTYPE:ENV  35.34972  6.574869  5.376490     P   0
## units!units  179.14713  7.110011 25.196462     P   0
## units!R        1.00000        NA        NA     F   0
## [1] "------- Analyzing trait DON -------"
## [1] "DON has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:26 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -21231.93           1.0   2196 13:35:26    0.0
##  2     -17367.32           1.0   2196 13:35:26    0.0
##  3     -13048.60           1.0   2196 13:35:26    0.0
##  4      -9918.73           1.0   2196 13:35:26    0.0
##  5      -8036.34           1.0   2196 13:35:26    0.0
##  6      -7256.48           1.0   2196 13:35:26    0.0
##  7      -6932.90           1.0   2196 13:35:26    0.0
##  8      -6810.85           1.0   2196 13:35:26    0.0
##  9      -6770.62           1.0   2196 13:35:26    0.0
## 10      -6760.85           1.0   2196 13:35:26    0.0
## 11      -6758.98           1.0   2196 13:35:26    0.0
## 12      -6758.57           1.0   2196 13:35:26    0.0
## 13      -6758.52           1.0   2196 13:35:26    0.0
## 14      -6758.52           1.0   2196 13:35:26    0.0
##              component std.error   z.ratio bound %ch
## ENV           77.14525 49.088534  1.571553     P 0.1
## GENOTYPE     102.31633 12.302329  8.316826     P 0.0
## GENOTYPE:ENV  52.62786  5.225355 10.071633     P 0.0
## units!units  102.48418  4.276448 23.964793     P 0.0
## units!R        1.00000        NA        NA     F 0.0
## [1] "------- Analyzing trait PH -------"
## [1] "PH has a single year of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May  7 13:35:26 2023
##           LogLik        Sigma2     DF     wall    cpu
##  1     -3067.202           1.0    842 13:35:26    0.0
##  2     -2656.660           1.0    842 13:35:26    0.0
##  3     -2215.257           1.0    842 13:35:26    0.0
##  4     -1928.060           1.0    842 13:35:26    0.0
##  5     -1793.587           1.0    842 13:35:26    0.0
##  6     -1765.862           1.0    842 13:35:26    0.0
##  7     -1763.549           1.0    842 13:35:26    0.0
##  8     -1763.517           1.0    842 13:35:26    0.0
##  9     -1763.517           1.0    842 13:35:26    0.0
##              component std.error    z.ratio bound %ch
## ENV           2.290131 3.3092394  0.6920414     P   0
## GENOTYPE     24.369899 3.0176516  8.0757828     P   0
## GENOTYPE:ENV  3.241033 1.0384778  3.1209455     P   0
## units!units  13.092888 0.8586791 15.2477070     P   0
## units!R       1.000000        NA         NA     F   0
kable(h2s, row.names = F, caption="Heritability Estimates and Their Resultant Standard Errors")
Heritability Estimates and Their Resultant Standard Errors
Trait Type Estimate SE
HD Per-Plot 0.5271801 0.0308343
HD Entry-Mean 0.9345059 0.0103703
FHB Per-Plot 0.5766447 0.0287133
FHB Entry-Mean 0.9570499 0.0066519
FDK Per-Plot 0.5680795 0.0280565
FDK Entry-Mean 0.9690813 0.0049837
DON Per-Plot 0.3974555 0.0304737
DON Entry-Mean 0.9018159 0.0132475
PH Per-Plot 0.5987128 0.0356496
PH Entry-Mean 0.8825275 0.0216826
write.csv(h2s,
          "heritabilities.csv",
          row.names = FALSE)

Make the rQTL file

#manipulate the genetic file to look like r qtl wants it
rownames<-colnames(geno)
geno<-as.data.frame(t(geno))
geno<-cbind(rownames,geno)
colnames(geno)<-geno[1,]
geno<-geno[-1,]
rownames(geno)<-c()
geno<-rownames_to_column(geno, var="dummy") #create a dummy column to sort later on
geno$dummy<-as.numeric(geno$dummy)
colnames(geno)[2]<-"GENOTYPE"

#merge all the blues
BLUEs_all<-left_join(BLUEs_mem, BLUEs_loc, by="GENOTYPE")

#merge the pheno and geno files
insertme<-full_join(BLUEs_all, geno, by="GENOTYPE")
insertme<-insertme[order(insertme$dummy),] #order by dummy
insertme<-insertme[1,]
pheno_geno<-merge(BLUEs_all, geno, by="GENOTYPE", all.x = T)
pheno_geno<-rbind(insertme,pheno_geno)
pheno_geno<-pheno_geno[order(pheno_geno$dummy),] #order by dummy
pheno_geno<-pheno_geno %>% filter(!is.na(dummy)==T) #get rid of any genotype that is not present
pheno_geno[1,(1:30)]<-"" #replace any NAs with blank data
dropme<-c("13955-AGS-2026","13955-JAMESTOWN", "13955-NCAG11", "13955-NC13-20076", "13955-GA06493-13LE6")
pheno_geno<-pheno_geno[!(pheno_geno$GENOTYPE %in% dropme),]
pheno_geno<-pheno_geno %>%
  select(-dummy)
colnames(pheno_geno)<-gsub("FHB", "VR", colnames(pheno_geno))

#write out rqtl file
write.csv(pheno_geno,
          "rqtl_input_file.csv",
          row.names = F)

Make a pairs plot of traits

pairs_plot_dat<-BLUEs_mem
colnames(pairs_plot_dat)<-c("Genotype",
                            "Visual Rating (1-9)",
                            "Fusarium Damaged Kernels (%)",
                            "Deoxynivalenol Content (PPM)",
                            "Heading Date (Days)",
                            "Plant Height (cm)")

jpeg(filename = "pairs_plot.jpg",
     width = 9,
     height = 9,
     units = "in",
     res = 320)

pairs.panels(pairs_plot_dat[,2:ncol(pairs_plot_dat)],
             hist.col = "gray",
             lm = TRUE,
             stars = TRUE,
             digits = 2,
             density = FALSE,
             ellipses = FALSE)

dev.off()
## png 
##   2

Build Linkage Map

#read in file
cross_file<-read.cross("csv", 
                       file="rqtl_input_file.csv",
                       genotypes=c("A","B","-"), 
                       alleles=c("A","B"), 
                       crosstype="dh")
##  --Read the following data:
##   187  individuals
##   2954  markers
##   29  phenotypes
## Warning in summary.cross(cross): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
##   (Perhaps it is in basepairs?)
##  --Cross type: dh
#Pull distorted markers
gt<-geno.table(cross_file)
dropme<-rownames(gt[gt$P.value<0.0001,])
cross_file<-drop.markers(cross_file, dropme)

#remove junk
remove(dropme, gt)
#create map with cM instead of BP
cross_file<- mstmap(cross_file, 
                    pop.type="dh", 
                    id = "GENOTYPE",
                    chr=c("1A", 
                          "1B", 
                          "1D", 
                          "2A", 
                          "2B", 
                          "2D", 
                          "3A", 
                          "3B", 
                          "3D", 
                          "4A", 
                          "4B", 
                          "4D", 
                          "5A", 
                          "5B", 
                          "5D", 
                          "6A", 
                          "6B", 
                          "6D", 
                          "7A", 
                          "7B", 
                          "7D"), 
                    anchor = TRUE, 
                    detectBadData = TRUE,
                    bychr = TRUE, 
                    miss.thresh = .15, 
                    mvest.bc = FALSE)
#jitter map to improve marker distance
cross_file<-jittermap(cross_file)
summary(cross_file)
##     Doubled haploids
## 
##     No. individuals:    187 
## 
##     No. phenotypes:     29 
##     Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100 
##                         100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9 
##                         100 100 100 100 
## 
##     No. chromosomes:    85 
##         Autosomes:      1A.1 1A.2 1B.1 1B.2 1D.1 1D.2 1D.3 1D.4 1D.5 1D.6 2A.1 
##                         2A.2 2A.3 2A.4 2A.5 2B.1 2B.2 2B.3 2B.4 2D.1 2D.2 2D.3 
##                         2D.4 2D.5 2D.6 2D.7 2D.8 3A.1 3A.2 3A.3 3A.4 3B.1 3B.2 
##                         3B.3 3B.4 3B.5 3D.1 3D.2 3D.3 3D.4 3D.5 4A.1 4A.2 4B.1 
##                         4B.2 4B.3 4B.4 4D.1 4D.2 5A.1 5A.2 5A.3 5A.4 5B.1 5B.2 
##                         5D.1 5D.2 5D.3 5D.4 6A.1 6A.2 6A.3 6B.1 6B.2 6B.3 6D.1 
##                         6D.2 7A.1 7A.2 7A.3 7A.4 7A.5 7B.1 7B.2 7B.3 7B.4 7B.5 
##                         7D.1 7D.2 7D.3 7D.4 7D.5 7D.6 7D.7 7D.8 
## 
##     Total markers:      2248 
##     No. markers:        111 1 133 1 57 1 1 1 1 1 133 1 1 1 1 137 4 1 1 27 1 2 1 
##                         1 2 1 4 177 1 1 1 218 1 1 1 1 1 14 1 1 1 120 4 70 1 1 1 
##                         1 12 154 1 1 1 97 5 18 1 1 4 124 2 1 123 1 1 50 2 221 1 
##                         1 1 1 123 1 1 1 1 40 1 3 1 1 2 3 1 
##     Percent genotyped:  92 
##     Genotypes (%):      AA:54.8  BB:45.2
names(cross_file$geno)
##  [1] "1A.1" "1A.2" "1B.1" "1B.2" "1D.1" "1D.2" "1D.3" "1D.4" "1D.5" "1D.6"
## [11] "2A.1" "2A.2" "2A.3" "2A.4" "2A.5" "2B.1" "2B.2" "2B.3" "2B.4" "2D.1"
## [21] "2D.2" "2D.3" "2D.4" "2D.5" "2D.6" "2D.7" "2D.8" "3A.1" "3A.2" "3A.3"
## [31] "3A.4" "3B.1" "3B.2" "3B.3" "3B.4" "3B.5" "3D.1" "3D.2" "3D.3" "3D.4"
## [41] "3D.5" "4A.1" "4A.2" "4B.1" "4B.2" "4B.3" "4B.4" "4D.1" "4D.2" "5A.1"
## [51] "5A.2" "5A.3" "5A.4" "5B.1" "5B.2" "5D.1" "5D.2" "5D.3" "5D.4" "6A.1"
## [61] "6A.2" "6A.3" "6B.1" "6B.2" "6B.3" "6D.1" "6D.2" "7A.1" "7A.2" "7A.3"
## [71] "7A.4" "7A.5" "7B.1" "7B.2" "7B.3" "7B.4" "7B.5" "7D.1" "7D.2" "7D.3"
## [81] "7D.4" "7D.5" "7D.6" "7D.7" "7D.8"
#remove small linkage groups of a specific length
chrLengths<-c() #set up list for chromosome lengths

for (i in cross_file$geno){
  chrLengths <- c(chrLengths, length(i$map))
}

#create data frame with number of markers per linkage group
chrLengths<-as.data.frame(cbind(names(cross_file$geno), chrLengths))
colnames(chrLengths)<-c("chr", "n_markers") 
chrLengths$n_markers<-as.numeric(chrLengths$n_markers) 
chrLengths<-filter(chrLengths, n_markers > 5) 

#new number of linkage groups
nrow(chrLengths)
## [1] 21
chrLengths
##     chr n_markers
## 1  1A.1       111
## 2  1B.1       133
## 3  1D.1        57
## 4  2A.1       133
## 5  2B.1       137
## 6  2D.1        27
## 7  3A.1       177
## 8  3B.1       218
## 9  3D.2        14
## 10 4A.1       120
## 11 4B.1        70
## 12 4D.2        12
## 13 5A.1       154
## 14 5B.1        97
## 15 5D.1        18
## 16 6A.1       124
## 17 6B.1       123
## 18 6D.1        50
## 19 7A.1       221
## 20 7B.1       123
## 21 7D.1        40
#subset the map with the linkage groups that pass the filter
cross_file<-subset(cross_file, chr = chrLengths$chr)
cross_file<-jittermap(cross_file) 

#show summary and write file with output
summary(cross_file)
##     Doubled haploids
## 
##     No. individuals:    187 
## 
##     No. phenotypes:     29 
##     Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100 
##                         100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9 
##                         100 100 100 100 
## 
##     No. chromosomes:    21 
##         Autosomes:      1A.1 1B.1 1D.1 2A.1 2B.1 2D.1 3A.1 3B.1 3D.2 4A.1 4B.1 
##                         4D.2 5A.1 5B.1 5D.1 6A.1 6B.1 6D.1 7A.1 7B.1 7D.1 
## 
##     Total markers:      2159 
##     No. markers:        111 133 57 133 137 27 177 218 14 120 70 12 154 97 18 124 
##                         123 50 221 123 40 
##     Percent genotyped:  92 
##     Genotypes (%):      AA:54.7  BB:45.3
write.cross(cross_file, "cross_file_without_parents",format = "csv")

#remove junk
remove(i, chrLengths)

#check map metrics
names(cross_file$pheno)
##  [1] "GENOTYPE"    "VR_ME"       "FDK_ME"      "DON_ME"      "HD_ME"      
##  [6] "PH_ME"       "VR_KIN2019"  "FDK_KIN2019" "DON_KIN2019" "HD_KIN2020" 
## [11] "VR_KIN2020"  "FDK_KIN2020" "DON_KIN2020" "PH_KIN2020"  "HD_LW2019"  
## [16] "VR_LW2019"   "FDK_LW2019"  "DON_LW2019"  "HD_LW2020"   "VR_LW2020"  
## [21] "FDK_LW2020"  "DON_LW2020"  "PH_LW2020"   "FDK_VIR2019" "DON_VIR2019"
## [26] "HD_VIR2020"  "VR_VIR2020"  "FDK_VIR2020" "DON_VIR2020"
names(cross_file)
## [1] "geno"  "pheno"
names(cross_file$geno)
##  [1] "1A.1" "1B.1" "1D.1" "2A.1" "2B.1" "2D.1" "3A.1" "3B.1" "3D.2" "4A.1"
## [11] "4B.1" "4D.2" "5A.1" "5B.1" "5D.1" "6A.1" "6B.1" "6D.1" "7A.1" "7B.1"
## [21] "7D.1"
nphe(cross_file)
## [1] 29
nind(cross_file)
## [1] 187
nchr(cross_file)
## [1] 21
summary(cross_file)
##     Doubled haploids
## 
##     No. individuals:    187 
## 
##     No. phenotypes:     29 
##     Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100 
##                         100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9 
##                         100 100 100 100 
## 
##     No. chromosomes:    21 
##         Autosomes:      1A.1 1B.1 1D.1 2A.1 2B.1 2D.1 3A.1 3B.1 3D.2 4A.1 4B.1 
##                         4D.2 5A.1 5B.1 5D.1 6A.1 6B.1 6D.1 7A.1 7B.1 7D.1 
## 
##     Total markers:      2159 
##     No. markers:        111 133 57 133 137 27 177 218 14 120 70 12 154 97 18 124 
##                         123 50 221 123 40 
##     Percent genotyped:  92 
##     Genotypes (%):      AA:54.7  BB:45.3
#rename linkage groups
names(cross_file$geno)<-c("1A", 
                          "1B", 
                          "1D", 
                          "2A", 
                          "2B", 
                          "2D", 
                          "3A", 
                          "3B", 
                          "3D", 
                          "4A", 
                          "4B", 
                          "4D",
                          "5A", 
                          "5B", 
                          "5D",
                          "6A", 
                          "6B", 
                          "6D", 
                          "7A", 
                          "7B", 
                          "7D")

#plot images of map
plotMap(cross_file, horizontal = FALSE, shift = FALSE, main="NC13-20076 Double Haploid Linkage Map")

geno.image(cross_file, main = "")

Check Map

#check orentation of markers
mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)
## $Markers
## [1] "character"
## 
## $Chr
## [1] "character"
## 
## $BP
## [1] "character"
## 
## $cM
## [1] "numeric"
pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
  geom_point()+
  facet_wrap(~Chr,
             ncol = 3,
             scales = "free")+
  labs(title = "Centimorgan Position vs. Megabase Pair",
       x = "Centimorgan (cM) Position",
       y = "Megabase Pair (Mbp) Position")

for (i in names(cross_file$geno)){
  q<-subset(pos, pos$Chr==i)
  if (q[1,3]>q[(nrow(q)),3]){
    cross_file<-flip.order(cross_file, i)
    print(paste("Chromsome",i,"has been flipped"))
    remove(q)
  }else{
    print(paste("Chromsome",i,"is in correct order"))
    remove(q)
  }
}
## [1] "Chromsome 1A is in correct order"
## [1] "Chromsome 1B is in correct order"
## [1] "Chromsome 1D is in correct order"
## [1] "Chromsome 2A is in correct order"
## [1] "Chromsome 2B is in correct order"
## [1] "Chromsome 2D is in correct order"
## [1] "Chromsome 3A is in correct order"
## [1] "Chromsome 3B is in correct order"
## [1] "Chromsome 3D is in correct order"
## [1] "Chromsome 4A is in correct order"
## [1] "Chromsome 4B is in correct order"
## [1] "Chromsome 4D is in correct order"
## [1] "Chromsome 5A is in correct order"
## [1] "Chromsome 5B is in correct order"
## [1] "Chromsome 5D is in correct order"
## [1] "Chromsome 6A is in correct order"
## [1] "Chromsome 6B is in correct order"
## [1] "Chromsome 6D is in correct order"
## [1] "Chromsome 7A is in correct order"
## [1] "Chromsome 7B is in correct order"
## [1] "Chromsome 7D has been flipped"
mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)
## $Markers
## [1] "character"
## 
## $Chr
## [1] "character"
## 
## $BP
## [1] "character"
## 
## $cM
## [1] "numeric"
pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
  geom_point()+
  facet_wrap(~Chr,
             ncol = 3,
             scales = "free")+
  labs(title = "Centimorgan Position vs. Megabase Pair",
       x = "Centimorgan (cM) Position",
       y = "Megabase Pair (Mbp) Position")

#remove bad markers
removeme<-c("S1A_527395239", "S1A_208718176", "S7A_621573544")
cross_file<-drop.markers(cross_file, removeme)

#check one last time
mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)
## $Markers
## [1] "character"
## 
## $Chr
## [1] "character"
## 
## $BP
## [1] "character"
## 
## $cM
## [1] "numeric"
pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
  geom_point()+
  facet_wrap(~Chr,
             ncol = 3,
             scales = "free")+
  labs(title = "Centimorgan Position vs. Megabase Pair",
       x = "Centimorgan (cM) Position",
       y = "Megabase Pair (Mbp) Position")

ggsave(plot = last_plot(),
       filename = "cM_vs_bp_comparison.jpg",
       units = "in",
       width = 8,
       height = 12,
       dpi = 320)

#remove
remove(mn, pos, removeme)
#calculate marker probabilities
cross_file<-calc.genoprob(cross_file, 
                          step=2.0, 
                          off.end=0.0, 
                          error.prob=1.0e-4, 
                          map.function="kosambi",
                          stepwidth="fixed")

#calculate marker probabilities for a simulated genotype
cross_file<-sim.geno(cross_file, 
                     n.draws = 128, 
                     step = 2, 
                     off.end = 0.0, 
                     error.prob=1.0e-4, 
                     map.function="kosambi", 
                     stepwidth="fixed")

#detect the number of cores of the computer processor 
ncores<-as.numeric(detectCores())

#set the number of permutations
nperms=1000

#### Run Heading Date ####
traits<-names(cross_file$pheno)[c(grep("HD", names(cross_file$pheno)),grep("PH", names(cross_file$pheno)))]

#make results vector
results<-list()

#run initial interval mapping and pull out QTL
for (i in traits){
  
  #announce 
  print(paste("------------ Interval Mapping of", i,"------------"))
  
  
  #perform IM with multiple imputation method
  print("Interval mapping...")
  scans<-scanone(cross_file, 
                 pheno.col = i,
                 addcovar = NULL,
                 model = "normal", 
                 method = "hk") 
  print("Done")
  
  #perform IM permutations mapping to define significance threshold
  print("Permutational interval mapping...")
  perms<-scanone(cross_file, 
                 pheno.col = i, 
                 addcovar = NULL,
                 model = "normal", 
                 method = "hk", 
                 n.perm = nperms, 
                 n.cluster = ncores-1) #set threshold
  print("Done")
  
  #plot QTL Scan
  print("Plotting...")
  threshold<-summary(perms, alpha=0.05)
  plot(scans,main=paste("IM for", i)) 
  abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
  legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
  
  #print plot
  jpeg(paste("Scan_IM_", i, ".jpg", sep = ""),
       width = 11,
       height = 4,
       units = "in",
       res = 320)
  threshold<-summary(perms, alpha=0.05)
  plot(scans,main=paste("IM for", i)) 
  abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
  legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
  dev.off()
  print("Done")
  
  
  #show the peak markers for QTL
  print("Defining QTL...")
  qtl<-summary(scans, perm=perms, 
               lodcolum=1, alpha=0.05)
  
  #rename QTL to identify which scan they came from
  qtl$name=paste("IM_", qtl$chr,"-pos-",qtl$pos,sep="")
  print("Done")
  
  #Define the QTL locations and effects
  print("Drawing QTL...")
  colnames(scans)<-c("chr","pos","lod")
  
  #place outputs in lists
  results$IM$scan[[i]]<-scans #place the scan a the list
  results$IM$perms[[i]]<-perms #place the permutations a the list
  results$IM$threshold[[i]]<-threshold #place the thresholds a the list
  
  #set up objects for defining QTL
  c<-qtl[,1] #define the chromosomes where QTL are found
  p<-qtl[,2] #define the positions of the QTL
  a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
  
  #Make new cross with genome subset
  a<-sim.geno(a, 
              n.draws = 128, 
              step = 2, 
              off.end = 0.0, 
              error.prob=1.0e-4, 
              map.function="kosambi", 
              stepwidth="fixed")
  
  #make a QTL object from that subset
  madeqtl<-makeqtl(a, 
                   c, 
                   p, 
                   qtl.name = qtl[,4], 
                   what = c("prob")) 

  #place that QTL object in a list
  results$IM$qtl[[i]]<-madeqtl 
  print("Done")
  
  #announce
  print("Running drop-one QTL analysis to check significance...")
  
  #test the significance of those QTL using drop one analysis
  qtlfit<-fitqtl(cross_file, 
                 qtl=results$IM$qtl[[i]], 
                 pheno.col = i,
                 model="normal", 
                 method = "hk") 
  
  #remove insignificant qtl
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))
    madeqtl<-dropfromqtl(madeqtl,
                         qtl.name = a$name)
    
    results$IM$qtl[[i]]<-madeqtl 
    
    remove(a,c,p,scans,perms,threshold)
    
  }else{
    
    remove(a,c,p,scans,perms,threshold)
    
  }
  
  #set vectors outside the loop
  j=1
  
  #make initial check object
  qtl_check<-results$IM$qtl[[i]]
  
  #run MQM until no significant peaks!
  repeat{
    
    #announce  
    print(paste("-------- Performing Additional QTL Scan for Trait", i,"--------"))
    
    mqm<-addqtl(cross = cross_file, 
                pheno.col = i, 
                qtl=qtl_check, 
                covar = NULL, 
                method = "hk")
    results[[paste("MQM", j, sep="")]]$scan[[i]]<-mqm
    
    #plot QTL Scan
    print("Plotting...")
    plot(mqm,main=paste("MQM", j, "for", i)) 
    abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
    legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
    print("Done")
    
    #write out picture to directory
    jpeg(paste("Scan_",paste("MQM", j, sep=""),"_", i, ".jpg", sep = ""),
     width = 11,
     height = 4,
     units = "in",
     res = 320)
    plot(mqm,main=paste("MQM", j, "for", i)) 
    abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
    legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
    dev.off()

    #make a qtl object
    qtl<-summary(mqm, 
                 perm=results$IM$perms[[i]], 
                 lodcolum=1, 
                 alpha=0.05)
  
    if(nrow(qtl)==0){
      
      results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
      print(paste("There were no new QTL identified for", i, "aborting loop"))
      break
      
    }else{
      
      #rename QTL to identify which scan they came from
      qtl$name=paste(paste("MQM", j, "_", sep=""), qtl$chr,"-pos-",qtl$pos,sep="")
       
      #Defining QTL
      print("Drawing new QTL")
      c<-qtl[,1] #define the chromosomes where QTL are found
      p<-qtl[,2] #define the positions of the QTL
      a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
       
      #simulate the genome for that subset
      a<-sim.geno(a, 
                  n.draws = 128, 
                  step = 2, 
                  off.end = 0.0, 
                  error.prob=1.0e-4, 
                  map.function="kosambi", 
                  stepwidth="fixed") 
       
      #make a QTL object from that subset
      madeqtl<-makeqtl(a, 
                       c, 
                       p, 
                       qtl.name = qtl[,4], 
                       what = c("prob")) 
   
      #place that QTL object in a list
      results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
      print("Done")
      
      #announce 
      print("Running drop-one QTL analysis...")
       
      #test the significance of those QTL using dropone analysis
      qtlfit<-fitqtl(cross_file, 
                     qtl=madeqtl, 
                     pheno.col = i, 
                     get.ests = T,
                     model="normal", 
                     method = "hk") 
      
      if(!is.null(qtlfit$result.drop)){
      
        print("Checking drop-one analysis for significance")
        a<-as.data.frame(qtlfit$result.drop)
        a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
        a<-data.frame(name=rownames(a[a$sig==0,]))
        a<-tidyr::separate(data = a, 
                           col = "name",
                           into = c("chr", "trash", "pos"),
                           sep = "-")
        a$pos<-as.numeric(a$pos)
        a$chr<-gsub(paste("MQM", j, "_", sep=""),"",a$chr)
        madeqtl<-dropfromqtl(madeqtl,
                             chr = a$chr,
                             pos = a$pos)
        
        madeqtl<-addtoqtl(cross_file,
                          qtl = madeqtl,
                          chr = qtl_check$chr,
                          pos = qtl_check$pos,
                          qtl.name = qtl_check$name)
              
        
        results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
        
        remove(a,c,p)
        
        
        }else{
          
          print(paste("There was only one new QTL identified for ",
                      i,
                      " in ", 
                      paste("MQM", j, sep=""),
                      "... Checking significance!",
                      sep = ""))
          
          if(qtlfit$result.full[1,6]<0.05){
      
            print("New QTL is significant, placing in total model!")
            
            madeqtl<-addtoqtl(cross_file,
                              qtl = madeqtl,
                              chr = qtl_check$chr,
                              pos = qtl_check$pos,
                              qtl.name = qtl_check$name)
            
            results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
      
            }else{
    
              results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL 
              print(paste("There were no new QTL identified for", i, "aborting loop!"))
              break
      
              }
          }
      }
      
    
    qtl_check<-results[[paste("MQM", j, sep="")]]$qtl[[i]]
    
    j=j+1
    
  }
  
  #check
  if(j-1==0){
    
    #pull and make final qtl object
    finalqtl<-results$IM$qtl[[i]]
    
  }else{
  
    #pull and make final qtl object
    finalqtl<-results[[paste("MQM", j-1, sep="")]]$qtl[[i]]
      
  }
  
  #check for significance
  qtlfit<-fitqtl(cross_file, 
                 qtl=finalqtl, 
                 pheno.col = i, 
                 get.ests = T,
                 model="normal", 
                 method = "hk") 
  
  #rename the qtl
  a<-finalqtl
  a<-data.frame(summary(a))
  a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
  a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
  finalqtl$name=a$name
  remove(a)
  
  #print
  print("Filtering insignificant markers...")
  
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))
    finalqtl<-dropfromqtl(finalqtl,
                          qtl.name = a$name)
    
    results$final$qtl_unrefine[[i]]<-finalqtl 
    
    remove(a, finalqtl)
    
  }else{
    
    remove(a, finalqtl)
    
  }
   
  #print
  print("Refining QTL position...")
  
  #refine
  results$final$qtl_refine[[i]]<-refineqtl(cross = cross_file, 
                                           pheno.col = i, 
                                           qtl = results$final$qtl_unrefine[[i]], 
                                           verbose = FALSE,
                                           method = "hk")
  #rename the qtl
  a<-results$final$qtl_refine[[i]]
  a<-data.frame(summary(a))
  a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
  a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
  results$final$qtl_refine[[i]]$name=a$name
  remove(a)
  
  #check for significance
  qtlfit<-fitqtl(cross_file, 
                 qtl=results$final$qtl_refine[[i]], 
                 pheno.col = i, 
                 get.ests = T,
                 model="normal", 
                 method = "hk") 
  
  finalqtl<-results$final$qtl_refine[[i]]
  
  #print
  print("Filtering insignificant markers...")
  
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))

    finalqtl<-dropfromqtl(finalqtl,
                          qtl.name = a$name)
    
    results$final$qtl_refine_filt[[i]]<-finalqtl 
    
    remove(a, finalqtl)
    
  }else{
    
    remove(a, finalqtl)
    
  }
  
  droplodint<-c()
  
  #create QTL 1.5 drop LOD intervals
  for(z in results$final$qtl_refine_filt[[i]]$name){
    
    a<-data.frame(name=z)
    a<-tidyr::separate(data = a,
                       col = "name",
                       into = c("scan", "other"),
                       sep = "_",
                       remove = FALSE)
    a<-tidyr::separate(data = a,
                       col = "other",
                       into = c("chr", "trash", "pos"),
                       sep = "-",
                       remove = TRUE)
    b<-results[[a$scan]]$scan[[i]]
    c<-find.pseudomarker(cross = cross_file,
                         chr = a$chr,
                         pos = as.numeric(a$pos))
    d<-lodint(results = b,
              chr = a$chr,
              qtl.index = c,
              drop = 1.5,
              lodcolumn = 1,
              expandtomarkers = TRUE)
    
    d<-data.frame(name=z, 
                  trait=i,
                  chr=unique(d$chr),
                  start_marker=rownames(d)[1],
                  start=d[1,2], 
                  mid_marker=rownames(d)[2],
                  mid_pos=d[2,2],
                  stop_marker=rownames(d)[3],
                  stop=d[3,2])
    
    droplodint<-rbind(droplodint, d)
    
    remove(a,b,c,d)
  }
  
  #place in results
  results$final$qtl_support_intervals[[i]]<-droplodint
  
  #remove
  remove(droplodint)
  
}
## [1] "------------ Interval Mapping of HD_ME ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_ME --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_ME --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for HD_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_KIN2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_KIN2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_KIN2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_KIN2020 --------"
## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for HD_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_LW2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = NULL, : Dropping 2 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "There was only one new QTL identified for HD_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = NULL, : Dropping 2 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for HD_LW2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.

## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_LW2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for HD_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_VIR2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_VIR2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_VIR2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_VIR2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for HD_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_ME ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_ME --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for PH_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait PH_ME --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for PH_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_KIN2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for PH_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_LW2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."

## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_LW2020 --------"
## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for PH_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
#rename results
results_HD_PH<-results

Pull Marker Covariates

covars<-rbind(results_HD_PH$final$qtl_support_intervals$HD_ME,
              results_HD_PH$final$qtl_support_intervals$PH_ME)

covar_markers<-data.frame(GENOTYPE=cross_file$pheno$GENOTYPE)

for(i in 1:nrow(covars)){
  
  a<-covars[i,]
  
  b<-fill.geno(cross = cross_file,
               method = c("no_dbl_XO"),
               map.function = "kosambi")
  
  b<-pull.geno(cross = b,
               chr = a$chr)
  
  a<-find.marker(cross = cross_file,
                 chr = a$chr,
                 pos = a$mid_pos)
  
  b<-as.data.frame(b[,a])
  
  colnames(b)=a
  
  covar_markers<-cbind(covar_markers, b)
  
  remove(a,b)
  
}

Visualize HD and PH QTL

#make object for ggplot
traits<-names(cross_file$pheno)[c(grep("HD", names(cross_file$pheno)),grep("PH", names(cross_file$pheno)))]

#make object
hd_ph_ggplot<-c()

for(i in traits){
  
  a<-results_HD_PH$final$qtl_support_intervals[[i]]
  hd_ph_ggplot<-rbind(hd_ph_ggplot, a)
  remove(a)
  
}

d<-results_HD_PH$IM$scan$HD_ME %>%
  filter(chr %in% hd_ph_ggplot$chr)

hd_ph_ggplot<-hd_ph_ggplot %>% 
  mutate(chr=paste(chr, trait, sep = "_"))

library(ggrepel)
library(patchwork)

finalplot<-list()
j=1

for(i in unique(d$chr)){

  a<-hd_ph_ggplot[grep(i, hd_ph_ggplot$chr),]
  a$trait<-factor(a$trait, levels = unique(hd_ph_ggplot$trait))
  b<-d[d$chr==i,]
  b<-b[-grep("loc", rownames(b)),]
  b<-rownames_to_column(b, var = "marker")
  Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
  
  if(j==1){
  
    c<-ggplot(data=b, aes(x=chr, y=pos))+
    geom_point(size=2)+
    geom_line(linewidth=1)+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("HD_ME" = "blue",
                                  "HD_KIN2020" = "dodgerblue",
                                  "HD_LW2019" = "darkslategray",
                                  "HD_LW2020" = "cyan",
                                  "HD_VIR2020" = 'cadetblue',
                                  "PH_ME" = 'violet',
                                  "PH_KIN2020" = 'darkorchid',
                                  "PH_LW2020"="mediumorchid"),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank())+
    labs(title = i,
         color = "Trait Within Environment",
         y = "cM")
    
    if(length(Label)>0){
      c<-c+
        geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
    }
  
  }else{
    
    c<-ggplot(data = b, aes(x=chr, y=pos))+
    geom_point(size=2)+
    geom_line(linewidth=1)+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("HD_ME" = "blue",
                                  "HD_KIN2020" = "dodgerblue",
                                  "HD_LW2019" = "darkslategray",
                                  "HD_LW2020" = "cyan",
                                  "HD_VIR2020" = 'cadetblue',
                                  "PH_ME" = 'violet',
                                  "PH_KIN2020" = 'darkorchid',
                                  "PH_LW2020"="mediumorchid"),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank(),
          axis.line.y = element_blank(),
          axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank())+
    labs(title = i,
         color = "Trait Within Environment",
         y = "cM")
    
    if(length(Label)>0){
      c<-c+
        geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
    }
  
    
  }
  finalplot[[paste("chr",i,sep="_")]]<-c
  remove(a,b,c)
  j=j+1
    
}

a<-(finalplot$chr_1B+finalplot$chr_2A+finalplot$chr_3A+
  finalplot$chr_4A+finalplot$chr_5A+finalplot$chr_6A+
  finalplot$chr_7A+finalplot$chr_7B+finalplot$chr_7D)+
  plot_layout(guides = "collect", nrow = 1)+
  plot_annotation(title = "1.5-LOD Support Intervals",
                  subtitle = "Plant Height (PH) and Heading Date (HD)")

a

ggsave(plot = a,
       filename = "support_intervals_hd_and_ph.jpg",
       width = 14,
       height = 8,
       dpi = 320,
       units = "in")

Visualize HD and PH QTL for ME only

#make object for ggplot
traits<-c("HD_ME", "PH_ME")

#make object
hd_ph_ggplot<-c()

for(i in traits){
  
  a<-results_HD_PH$final$qtl_support_intervals[[i]]
  hd_ph_ggplot<-rbind(hd_ph_ggplot, a)
  remove(a)
  
}

d<-results_HD_PH$IM$scan$HD_ME %>%
  filter(chr %in% hd_ph_ggplot$chr)

hd_ph_ggplot<-hd_ph_ggplot %>% 
  mutate(trait=ifelse(trait=="HD_ME", "Heading Date", 
                      ifelse(trait=="PH_ME", "Plant Height", "ERROR"))) %>%
  mutate(chr=paste(chr, trait, sep = "_"))

library(ggrepel)
library(patchwork)

finalplot<-list()
j=1

for(i in unique(d$chr)){

  a<-hd_ph_ggplot[grep(i, hd_ph_ggplot$chr),]
  a$trait<-factor(a$trait, levels = unique(hd_ph_ggplot$trait))
  b<-d[d$chr==i,]
  b<-b[-grep("loc", rownames(b)),]
  b<-rownames_to_column(b, var = "marker")
  Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
  
  if(j==1){
  
    c<-ggplot(data=b, aes(x=chr, y=pos))+
    geom_point(size=2)+
    geom_line(linewidth=1)+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("Heading Date" = "blue",
                                  "Plant Height" = 'violet'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank())+
    labs(title = i,
         color = "Trait",
         y = "cM")
    
    if(length(Label)>0){
      c<-c+
        geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
    }
  
  }else{
    
    c<-ggplot(data = b, aes(x=chr, y=pos))+
    geom_point(size=2)+
    geom_line(linewidth=1)+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("Heading Date" = "blue",
                                  "Plant Height" = 'violet'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank(),
          axis.line.y = element_blank(),
          axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank())+
    labs(title = i,
         color = "Trait",
         y = "cM")
    
    if(length(Label)>0){
      c<-c+
        geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
    }
  
    
  }
  finalplot[[paste("chr",i,sep="_")]]<-c
  remove(a,b,c)
  j=j+1
    
}

a<-(finalplot$chr_2A+finalplot$chr_4A+finalplot$chr_5A+
    finalplot$chr_6A+finalplot$chr_7B+finalplot$chr_7D)+
  plot_layout(guides = "collect", nrow = 1)+
  plot_annotation(title = "1.5-LOD Support Intervals",
                  subtitle = "Plant Height (PH) and Heading Date (HD)")

a

ggsave(plot = a,
       filename = "support_intervals_hd_and_ph_ME_only.jpg",
       width = 14,
       height = 8,
       dpi = 320,
       units = "in")

Run Mapping for Disease Traits

#set the number of permutations
nperms=1000

#### Run disease traits ####
traits<-names(cross_file$pheno)[c(grep("VR", names(cross_file$pheno)),grep("FDK", names(cross_file$pheno)),grep("DON", names(cross_file$pheno)))]

#make results vector
results<-list()

#run initial interval mapping and pull out QTL
for (i in traits){
  
  #announce 
  print(paste("------------ Interval Mapping of", i,"------------"))
  
  
  #perform IM with multiple imputation method
  print("Interval mapping...")
  scans<-scanone(cross_file, 
                 pheno.col = i,
                 addcovar = covar_markers[,2:ncol(covar_markers)],
                 model = "normal", 
                 method = "hk") 
  print("Done")
  
  #perform IM permutations mapping to define significance threshold
  print("Permutational interval mapping...")
  perms<-scanone(cross_file, 
                 pheno.col = i, 
                 addcovar = covar_markers[,2:ncol(covar_markers)],
                 model = "normal", 
                 method = "hk", 
                 n.perm = nperms, 
                 n.cluster = ncores-1) #set threshold
  print("Done")
  
  #plot QTL Scan
  print("Plotting...")
  threshold<-summary(perms, alpha=0.05)
  plot(scans,main=paste("IM for", i)) 
  abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
  legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
  
  #print plot
  jpeg(paste("Scan_IM_", i, ".jpg", sep = ""),
       width = 11,
       height = 4,
       units = "in",
       res = 320)
  threshold<-summary(perms, alpha=0.05)
  plot(scans,main=paste("IM for", i)) 
  abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
  legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
  dev.off()
  print("Done")
  
  
  #show the peak markers for QTL
  print("Defining QTL...")
  qtl<-summary(scans, perm=perms, 
               lodcolum=1, alpha=0.05)
  
  #rename QTL to identify which scan they came from
  qtl$name=paste("IM_", qtl$chr,"-pos-",qtl$pos,sep="")
  print("Done")
  
  #Define the QTL locations and effects
  print("Drawing QTL...")
  colnames(scans)<-c("chr","pos","lod")
  
  #place outputs in lists
  results$IM$scan[[i]]<-scans #place the scan a the list
  results$IM$perms[[i]]<-perms #place the permutations a the list
  results$IM$threshold[[i]]<-threshold #place the thresholds a the list
  
  #set up objects for defining QTL
  c<-qtl[,1] #define the chromosomes where QTL are found
  p<-qtl[,2] #define the positions of the QTL
  a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
  
  #Make new cross with genome subset
  a<-sim.geno(a, 
              n.draws = 128, 
              step = 2, 
              off.end = 0.0, 
              error.prob=1.0e-4, 
              map.function="kosambi", 
              stepwidth="fixed")
  
  #make a QTL object from that subset
  madeqtl<-makeqtl(a, 
                   c, 
                   p, 
                   qtl.name = qtl[,4], 
                   what = c("prob")) 

  #place that QTL object in a list
  results$IM$qtl[[i]]<-madeqtl 
  print("Done")
  
  #announce
  print("Running drop-one QTL analysis to check significance...")
  
  #test the significance of those QTL using drop one analysis
  qtlfit<-fitqtl(cross_file, 
                 qtl=results$IM$qtl[[i]], 
                 pheno.col = i,
                 model="normal", 
                 method = "hk") 
  
  #remove insignificant qtl
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))
    madeqtl<-dropfromqtl(madeqtl,
                         qtl.name = a$name)
    
    results$IM$qtl[[i]]<-madeqtl 
    
    remove(a,c,p,scans,perms,threshold)
    
  }else{
    
    remove(a,c,p,scans,perms,threshold)
    
  }
  
  #set vectors outside the loop
  j=1
  
  #make initial check object
  qtl_check<-results$IM$qtl[[i]]
  
  #run MQM until no significant peaks!
  repeat{
    
    #announce  
    print(paste("-------- Performing Additional QTL Scan for Trait", i,"--------"))
    
    mqm<-addqtl(cross = cross_file, 
                pheno.col = i, 
                qtl=qtl_check, 
                covar = covar_markers[,2:ncol(covar_markers)], 
                method = "hk")
    results[[paste("MQM", j, sep="")]]$scan[[i]]<-mqm
    
    #plot QTL Scan
    print("Plotting...")
    plot(mqm,main=paste("MQM", j, "for", i)) 
    abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
    legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
    print("Done")
    
    #write out picture to directory
    jpeg(paste("Scan_",paste("MQM", j, sep=""),"_", i, ".jpg", sep = ""),
     width = 11,
     height = 4,
     units = "in",
     res = 320)
    plot(mqm,main=paste("MQM", j, "for", i)) 
    abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
    legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
    dev.off()

    #make a qtl object
    qtl<-summary(mqm, 
                 perm=results$IM$perms[[i]], 
                 lodcolum=1, 
                 alpha=0.05)
  
    if(nrow(qtl)==0){
      
      results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
      print(paste("There were no new QTL identified for", i, "aborting loop"))
      break
      
    }else{
      
      #rename QTL to identify which scan they came from
      qtl$name=paste(paste("MQM", j, "_", sep=""), qtl$chr,"-pos-",qtl$pos,sep="")
       
      #Defining QTL
      print("Drawing new QTL")
      c<-qtl[,1] #define the chromosomes where QTL are found
      p<-qtl[,2] #define the positions of the QTL
      a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
       
      #simulate the genome for that subset
      a<-sim.geno(a, 
                  n.draws = 128, 
                  step = 2, 
                  off.end = 0.0, 
                  error.prob=1.0e-4, 
                  map.function="kosambi", 
                  stepwidth="fixed") 
       
      #make a QTL object from that subset
      madeqtl<-makeqtl(a, 
                       c, 
                       p, 
                       qtl.name = qtl[,4], 
                       what = c("prob")) 
   
      #place that QTL object in a list
      results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
      print("Done")
      
      #announce 
      print("Running drop-one QTL analysis...")
       
      #test the significance of those QTL using dropone analysis
      qtlfit<-fitqtl(cross_file, 
                     qtl=madeqtl, 
                     pheno.col = i, 
                     get.ests = T,
                     model="normal", 
                     method = "hk") 
      
      if(!is.null(qtlfit$result.drop)){
      
        print("Checking drop-one analysis for significance")
        a<-as.data.frame(qtlfit$result.drop)
        a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
        a<-data.frame(name=rownames(a[a$sig==0,]))
        a<-tidyr::separate(data = a, 
                           col = "name",
                           into = c("chr", "trash", "pos"),
                           sep = "-")
        a$pos<-as.numeric(a$pos)
        a$chr<-gsub(paste("MQM", j, "_", sep=""),"",a$chr)
        madeqtl<-dropfromqtl(madeqtl,
                             chr = a$chr,
                             pos = a$pos)
        
        madeqtl<-addtoqtl(cross_file,
                          qtl = madeqtl,
                          chr = qtl_check$chr,
                          pos = qtl_check$pos,
                          qtl.name = qtl_check$name)
              
        
        results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
        
        remove(a,c,p)
        
        
        }else{
          
          print(paste("There was only one new QTL identified for ",
                      i,
                      " in ", 
                      paste("MQM", j, sep=""),
                      "... Checking significance!",
                      sep = ""))
          
          if(qtlfit$result.full[1,6]<0.05){
      
            print("New QTL is significant, placing in total model!")
            
            madeqtl<-addtoqtl(cross_file,
                              qtl = madeqtl,
                              chr = qtl_check$chr,
                              pos = qtl_check$pos,
                              qtl.name = qtl_check$name)
            
            results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl 
      
            }else{
    
              results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL 
              print(paste("There were no new QTL identified for", i, "aborting loop!"))
              break
      
              }
          }
      }
      
    
    qtl_check<-results[[paste("MQM", j, sep="")]]$qtl[[i]]
    
    j=j+1
    
  }
  
  #check
  if(j-1==0){
    
    #pull and make final qtl object
    finalqtl<-results$IM$qtl[[i]]
    
  }else{
  
    #pull and make final qtl object
    finalqtl<-results[[paste("MQM", j-1, sep="")]]$qtl[[i]]
      
  }
  
  #check for significance
  qtlfit<-fitqtl(cross_file, 
                 qtl=finalqtl, 
                 pheno.col = i, 
                 get.ests = T,
                 model="normal", 
                 method = "hk") 
  
  #rename the qtl
  a<-finalqtl
  a<-data.frame(summary(a))
  a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
  a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
  finalqtl$name=a$name
  remove(a)
  
  #print
  print("Filtering insignificant markers...")
  
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))
    finalqtl<-dropfromqtl(finalqtl,
                          qtl.name = a$name)
    
    results$final$qtl_unrefine[[i]]<-finalqtl 
    
    remove(a, finalqtl)
    
  }else{
    
    remove(a, finalqtl)
    
  }
   
  #print
  print("Refining QTL position...")
  
  #refine
  results$final$qtl_refine[[i]]<-refineqtl(cross = cross_file, 
                                           pheno.col = i, 
                                           qtl = results$final$qtl_unrefine[[i]], 
                                           verbose = FALSE,
                                           method = "hk")
  #rename the qtl
  a<-results$final$qtl_refine[[i]]
  a<-data.frame(summary(a))
  a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
  a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
  results$final$qtl_refine[[i]]$name=a$name
  remove(a)
  
  #check for significance
  qtlfit<-fitqtl(cross_file, 
                 qtl=results$final$qtl_refine[[i]], 
                 pheno.col = i, 
                 get.ests = T,
                 model="normal", 
                 method = "hk") 
  
  finalqtl<-results$final$qtl_refine[[i]]
  
  #print
  print("Filtering insignificant markers...")
  
  if(!is.null(qtlfit$result.drop)){
    
    a<-as.data.frame(qtlfit$result.drop)
    a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
    a<-data.frame(name=rownames(a[a$sig==0,]))

    finalqtl<-dropfromqtl(finalqtl,
                          qtl.name = a$name)
    
    results$final$qtl_refine_filt[[i]]<-finalqtl 
    
    remove(a, finalqtl)
    
  }else{
    
    remove(a, finalqtl)
    
  }
  
  droplodint<-c()
  
  #create QTL 1.5 drop LOD intervals
  for(z in results$final$qtl_refine_filt[[i]]$name){
    
    a<-data.frame(name=z)
    a<-tidyr::separate(data = a,
                       col = "name",
                       into = c("scan", "other"),
                       sep = "_",
                       remove = FALSE)
    a<-tidyr::separate(data = a,
                       col = "other",
                       into = c("chr", "trash", "pos"),
                       sep = "-",
                       remove = TRUE)
    b<-results[[a$scan]]$scan[[i]]
    c<-find.pseudomarker(cross = cross_file,
                         chr = a$chr,
                         pos = as.numeric(a$pos))
    d<-lodint(results = b,
              chr = a$chr,
              qtl.index = c,
              drop = 1.5,
              lodcolumn = 1,
              expandtomarkers = TRUE)
    
    d<-data.frame(name=z, 
                  trait=i,
                  chr=unique(d$chr),
                  start_marker=rownames(d)[1],
                  start=d[1,2], 
                  mid_marker=rownames(d)[2],
                  mid_pos=d[2,2],
                  stop_marker=rownames(d)[3],
                  stop=d[3,2])
    
    droplodint<-rbind(droplodint, d)
    
    remove(a,b,c,d)
  }
  
  #place in results
  results$final$qtl_support_intervals[[i]]<-droplodint
  
  #remove
  remove(droplodint)
  
}
## [1] "------------ Interval Mapping of VR_ME ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait VR_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_KIN2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait VR_KIN2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_KIN2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_KIN2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_KIN2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_LW2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "There was only one new QTL identified for VR_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_LW2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_LW2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for VR_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for VR_LW2020 in MQM2... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## Warning: Expected 3 pieces. Additional pieces discarded in 1 rows [1].
## [1] "------------ Interval Mapping of VR_VIR2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_VIR2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for VR_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_ME ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_KIN2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_KIN2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_KIN2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_LW2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "There was only one new QTL identified for FDK_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_LW2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_LW2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for FDK_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_VIR2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_VIR2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_VIR2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for FDK_VIR2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for FDK_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_ME ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for DON_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_ME --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_KIN2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait DON_KIN2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_KIN2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_KIN2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_KIN2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_LW2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_LW2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_LW2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for DON_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_VIR2019 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "There was only one new QTL identified for DON_VIR2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2019 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"
## [1] "There were no new QTL identified for DON_VIR2019 aborting loop"
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_VIR2020 ------------"
## [1] "Interval mapping..."
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.

## [1] "Done"
## [1] "Permutational interval mapping..."
##  -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2020 --------"
## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.

## [1] "Plotting..."
## [1] "Done"

## [1] "There were no new QTL identified for DON_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
#rename results
results_FHB_FDK_DON<-results

Visualize FHB, FDK, and DON QTL

#make object for ggplot
traits<-names(cross_file$pheno)[c(grep("VR", names(cross_file$pheno)),grep("FDK", names(cross_file$pheno)),grep("DON", names(cross_file$pheno)))]

#make object
fhb_fdk_don_ggplot<-c()

for(i in traits){
  
  a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
  fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
  remove(a)
  
}

d<-results_FHB_FDK_DON$IM$scan$VR_ME %>%
  filter(chr %in% as.character(unique(fhb_fdk_don_ggplot$chr)))

fhb_fdk_don_ggplot<-fhb_fdk_don_ggplot %>% 
  mutate(trait=gsub("FHB", "VR", trait)) %>%
  mutate(chr=paste(chr, trait, sep = "_"))

library(ggrepel)
library(patchwork)

finalplot<-list()
j=1

for(i in unique(d$chr)){

  a<-fhb_fdk_don_ggplot[grep(i, fhb_fdk_don_ggplot$chr),]
  a$trait<-factor(a$trait, levels = unique(fhb_fdk_don_ggplot$trait))
  b<-d[d$chr==i,]
  b<-b[-grep("loc", rownames(b)),]
  b<-rownames_to_column(b, var = "marker")
  Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
  
  if(j==1 | j==11){
  
    c<-ggplot(data=b, aes(x=chr, y=pos))+
    geom_point(size=1)+
    geom_line()+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("VR_ME" = "red",
                                  "VR_KIN2019" = "firebrick",
                                  "VR_KIN2020" = "indianred",
                                  "VR_LW2019" = "brown",
                                  "VR_LW2020" = 'orangered',
                                  "VR_VIR2020" = 'tomato',
                                  "FDK_ME" = 'yellow',
                                  "FDK_KIN2019" = 'darkgoldenrod',
                                  "FDK_KIN2020" = 'gold',
                                  "FDK_LW2019" = 'goldenrod',
                                  "FDK_LW2020" = 'lightgoldenrod',
                                  "FDK_VIR2019" = 'darkgoldenrod2',
                                  "FDK_VIR2020" = 'khaki4',
                                  "DON_ME" = 'green',
                                  "DON_KIN2019" = 'chartreuse3',
                                  "DON_KIN2020" = 'darkolivegreen',
                                  "DON_LW2019" = 'darkseagreen',
                                  "DON_LW2020" = 'olivedrab',
                                  "DON_VIR2019" = 'palegreen',
                                  "DON_VIR2020" = 'springgreen1'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank())+
    labs(title = i,
         color = "Trait Within Environment",
         y = "cM")
  
  }else{
    
    c<-ggplot(data = b, aes(x=chr, y=pos))+
    geom_point(size=1)+
    geom_line()+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("VR_ME" = "red",
                                  "VR_KIN2019" = "firebrick",
                                  "VR_KIN2020" = "indianred",
                                  "VR_LW2019" = "brown",
                                  "VR_LW2020" = 'orangered',
                                  "VR_VIR2020" = 'tomato',
                                  "FDK_ME" = 'yellow',
                                  "FDK_KIN2019" = 'darkgoldenrod',
                                  "FDK_KIN2020" = 'gold',
                                  "FDK_LW2019" = 'goldenrod',
                                  "FDK_LW2020" = 'lightgoldenrod',
                                  "FDK_VIR2019" = 'darkgoldenrod2',
                                  "FDK_VIR2020" = 'khaki4',
                                  "DON_ME" = 'green',
                                  "DON_KIN2019" = 'chartreuse3',
                                  "DON_KIN2020" = 'darkolivegreen',
                                  "DON_LW2019" = 'darkseagreen',
                                  "DON_LW2020" = 'olivedrab',
                                  "DON_VIR2019" = 'palegreen',
                                  "DON_VIR2020" = 'springgreen1'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank(),
          axis.line.y = element_blank(),
          axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank())+
    labs(title = i,
         color = "Trait Within Environment",
         y = "cM")
    
  }
  finalplot[[paste("chr",i,sep="_")]]<-c
  remove(a,b,c)
  j=j+1
    
}

a<-(finalplot$chr_1A+finalplot$chr_1B+finalplot$chr_1D+
      finalplot$chr_2A+finalplot$chr_2B+finalplot$chr_2D+
      finalplot$chr_3A+finalplot$chr_3B+finalplot$chr_3D+
      finalplot$chr_4A+finalplot$chr_4B+finalplot$chr_4D+
      finalplot$chr_5A+finalplot$chr_5B+finalplot$chr_5D+
      finalplot$chr_6A+finalplot$chr_6B+
      finalplot$chr_7A+finalplot$chr_7D)+
  plot_layout(guides = "collect", 
              nrow = 2,
              ncol = 10)+
  plot_annotation(title = "1.5-LOD Support Intervals",
                  subtitle = "Visual Ratings (VR), Fusarium Damaged Kernels (FDK), and Deoxynivalenol Content (DON)")

a

ggsave(plot = a,
       filename = "support_intervals_vr_fdk_and_don.jpg",
       width = 20,
       height = 8,
       dpi = 320,
       units = "in")

Redraw with just ME

#make object for ggplot
traits<-c("VR_ME", "FDK_ME", "DON_ME")

#make object
fhb_fdk_don_ggplot<-c()

for(i in traits){
  
  a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
  fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
  remove(a)
  
}

d<-results_FHB_FDK_DON$IM$scan$VR_ME %>%
  filter(chr %in% as.character(unique(fhb_fdk_don_ggplot$chr)))

fhb_fdk_don_ggplot<-fhb_fdk_don_ggplot %>% 
  mutate(trait=ifelse(trait=="VR_ME", "Visual Rating",
                      ifelse(trait=="FDK_ME", "Fusarium Damaged Kernels",
                             ifelse(trait=="DON_ME", "Deoxynivalenol Content", "ERROR")))) %>%
  mutate(chr=paste(chr, trait, sep = "_"))

library(ggrepel)
library(patchwork)

finalplot<-list()
j=1

for(i in c("2A", "3B", "4A", "4B", "6A", "7A")){

  a<-fhb_fdk_don_ggplot[grep(i, fhb_fdk_don_ggplot$chr),]
  a$trait<-factor(a$trait, levels = unique(fhb_fdk_don_ggplot$trait))
  b<-d[d$chr==i,]
  b<-b[-grep("loc", rownames(b)),]
  b<-rownames_to_column(b, var = "marker")
  Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
  
  if(j==1){
  
    c<-ggplot(data=b, aes(x=chr, y=pos))+
    geom_point(size=1)+
    geom_line()+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("Visual Rating" = "red",
                                  "Fusarium Damaged Kernels" = 'yellow',
                                  "Deoxynivalenol Content" = 'green'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank())+
    labs(title = i,
         color = "Trait",
         y = "cM")
  
  }else{
    
    c<-ggplot(data = b, aes(x=chr, y=pos))+
    geom_point(size=1)+
    geom_line()+
    scale_y_reverse(limits = c(max(d$pos),0))+
    scale_color_manual(values = c("Visual Rating" = "red",
                                  "Fusarium Damaged Kernels" = 'yellow',
                                  "Deoxynivalenol Content" = 'green'),
                       drop=FALSE)+
    geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
    theme_classic()+
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.line.x = element_blank(),
          axis.line.y = element_blank(),
          axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank())+
    labs(title = i,
         color = "Trait",
         y = "cM")
    
  }
  finalplot[[paste("chr",i,sep="_")]]<-c
  remove(a,b,c)
  j=j+1
    
}

a<-(finalplot$chr_2A+finalplot$chr_3B+finalplot$chr_4A+finalplot$chr_4B+finalplot$chr_6A+finalplot$chr_7A)+
  plot_layout(guides = "collect", 
              nrow = 1)+
  plot_annotation(title = "1.5-LOD Support Intervals",
                  subtitle = "Visual Ratings (VR), Fusarium Damaged Kernels (FDK), and Deoxynivalenol Content (DON)")

a

ggsave(plot = a,
       filename = "support_intervals_vr_fdk_and_don_ME_only.jpg",
       width = 10,
       height = 8,
       dpi = 320,
       units = "in")

Percent variance in ME

fhb_fdk_don_ggplot<-c()

for(i in traits){
  
  a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
  fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
  remove(a)
  
}

qtl_info<-c()

for(i in c("2A", "3B", "4A", "4B", "6A", "7A")){

  a<-fhb_fdk_don_ggplot[fhb_fdk_don_ggplot$chr==i,]
  
  vr_effects<-fitqtl(cross = cross_file,
                     pheno.col = "VR_ME",
                     covar = covar_markers[,2:ncol(covar_markers)],
                     qtl = results_FHB_FDK_DON$final$qtl_refine_filt$VR_ME,
                     get.ests = TRUE)
      
  fdk_effects<-fitqtl(cross = cross_file,
                      pheno.col = "FDK_ME",
                      covar = covar_markers[,2:ncol(covar_markers)],
                      qtl = results_FHB_FDK_DON$final$qtl_refine_filt$FDK_ME,
                      get.ests = TRUE)
        
  don_effects<-fitqtl(cross = cross_file,
                      pheno.col = "DON_ME",
                      covar = covar_markers[,2:ncol(covar_markers)],
                      qtl = results_FHB_FDK_DON$final$qtl_refine_filt$DON_ME,
                      get.ests = TRUE)
  
  b<-as.data.frame(vr_effects$result.drop)[rownames(vr_effects$result.drop) %in% a$name,]
  b<-rownames_to_column(b, var = "name")
  b<-b %>% left_join(a, by="name") %>% filter(trait=="VR_ME")
  
  c<-as.data.frame(fdk_effects$result.drop)[rownames(fdk_effects$result.drop) %in% a$name,]
  c<-rownames_to_column(c, var = "name")
  c<-c %>% left_join(a, by="name") %>% filter(trait=="FDK_ME")
  
  d<-as.data.frame(don_effects$result.drop)[rownames(don_effects$result.drop) %in% a$name,]
  d<-rownames_to_column(d, var = "name")
  d<-d %>% left_join(a, by="name") %>% filter(trait=="DON_ME")
  
  e<-as.data.frame(summary(vr_effects)$est)[rownames(summary(vr_effects)$est) %in% a$name,]
  e<-rownames_to_column(e, var = "name")
  b<-b %>% left_join(e, by="name") %>% filter(trait=="VR_ME")
  
  f<-as.data.frame(summary(fdk_effects)$est)[rownames(summary(fdk_effects)$est) %in% a$name,]
  f<-rownames_to_column(f, var = "name")
  c<-c %>% left_join(f, by="name") %>% filter(trait=="FDK_ME")
  
  g<-as.data.frame(summary(don_effects)$est)[rownames(summary(don_effects)$est) %in% a$name,]
  g<-rownames_to_column(g, var = "name")
  d<-d %>% left_join(g, by="name") %>% filter(trait=="DON_ME")
  
  h<-rbind(b,c,d)
  
  qtl_info<-rbind(qtl_info, h)

  remove(a,b,c,d,e,f,g,h,vr_effects,fdk_effects,don_effects)
}
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".
## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.
qtl_info<-qtl_info %>%
  arrange(chr, mid_pos, trait) %>%
  separate(start_marker, into=c("junk1", "start_bp"), sep = "_", remove = FALSE) %>%
  separate(stop_marker, into=c("junk2", "stop_bp"), sep = "_", remove = FALSE) %>%
  mutate(`Proximal Mbp Position`= round(as.numeric(start_bp)/1000000, 2),
         `Distal Mbp Position` = round(as.numeric(stop_bp)/1000000, 2),
         newname = paste("QFhb.nc-", chr, sep = ""),
         `%var`=`%var`/100) %>%
  select(newname, 
         trait, 
         chr, 
         start_marker, 
         start, 
         `Proximal Mbp Position`,
         stop_marker, 
         stop, 
         `Distal Mbp Position`,
         est, 
         SE, 
         `%var`) %>%
  rename(QTL=newname,
         Trait=trait,
         Chromosome=chr,
         `Proximal Flanking Marker`=start_marker,
         `Proximal cM Position`=start,
         `Distal Flanking Marker`=stop_marker,
         `Distal cM Position`= stop,
         `Estimated Effect`=est,
         `Standard Error`=SE,
         `Percent Variance`=`%var`)

write.csv(qtl_info,
          "final_qtl_information.csv",
          row.names = FALSE)